(CC-BY-NC-SA)

last edit 2015-09-13

Introduction

Working with spatial data in R I find myself quite often in the need to quickly visually check whether a certain analysis has produced reasonable results. There are two ways I usually do this. Either I:

  1. (sp)plot the data in R and then toggle back and forth between the static plots (I use RStudio) or
  2. save the data to the disk and then open in QGIS or similar to interactively examine the results.

Both these approaches are semi-optimal. Where option 1. is fine for a quick glance at a coarse patterns, it lacks the possibility to have a closer look into the results via zooming and panning. While option 2. provides the interactivity, the detour via the hard disk is annoying (at best), especially when fine-tuning and checking regularly.

I attended this years useR2015! conference in Aalborg (which was marvelous!) and attended the session on interactive graphics in R where Joe Cheng from RStudio presented the leaflet package. leaflet is great as it gives you a great deal of control over the individual map components. This, however, also means that in order to get some useful visualization of spatial objects, we need to do a fair bit of coding for the map components to show what we want (e.g. we would need to provide popup-text manually for each object that we want to map). What a GIS-like functionality would need is some default behaviour for different objects from the spatial universe.

This got me thinking and sparked my enthusiasm to write some wrapper functions for leaflet to provide at least very basic GIS-like interactive graphing capabilities that are directly accessible within RStudio (or the web browser, if you’re not using RStudio).

The result is a package called mapview. The main workhorse function is mapView() (mind the capital ‘V’) and is currently defined for:

General design

A call to mapView() will return an object of class mapview. This class has 2 slots

  • @object - a list of the objects that are displayed on the map. This means that this slot will contain the re-projected (and in the case of Raster objects possibly re-sampled) objects which enables tracing of the modifications that took place.
  • @map - the leaflet map. This is an S3 class object (see leaflet documentation for details on the specifics).

By default mapView() provides two base layers between which one can toggle (a standard OpenStreeMap layer and ESRI’s WorldImagery layer). Depending on the object and/or argument settings one or several layers are created (each with or without it’s own legend) that can also be toggled. Note that in order to render properly, all layers need to be re-projected to leaflet’s underlying web mercator projection

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs

which in the case of large Raster* objects or Spatial objects with lots of features can be time consuming.

In the following I would like to present a few use case scenarios that highlight the current capabilities of mapview (largely taken from the help files):

Raster data objects

NOTE: similar to spplot() for Raster* objects mapView() has a maxpixels = argument to avoid long rendering times for large rasters. This is currently set to a default value of 500000 (which produces acceptable rendering times on my machine) and a suitable warning is printed for larger Raster objects.

RasterLayer

## to install mapview use 
# library(devtools)
# install_github("environmentalinformatics-marburg/mapview")
library(mapview)
library(raster)

data(meuse.grid)
coordinates(meuse.grid) = ~x+y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
gridded(meuse.grid) = TRUE

## convert to RasterStack
meuse_rst <- stack(meuse.grid)

mapView(meuse_rst[[4]])

Actually, in the original meuse.grid data frame the 4th column is a factor with 3 levels. So when we convert the RasterLayer to a factor, we will see that we get a proper factorial layer with a suitable legend.

meuse_rst[[4]] <- as.factor(meuse_rst[[4]])
mapView(meuse_rst[[4]])

RasterStack/Brick

If you pass a RasterStack/Brick to mapView() this will create one map layer + legend for each RasterLayer. This is likely to cause problems for larger Stacks as there is currently no way to tell mapView() to only show the legend of the highlighted layer. This means that some of the legends won’t show in the viewer.

mapView(meuse_rst)

There is (to my knowledge) currently no way to provide this functionality in the leaflet package, though I have seen solutions for this in JavaScript which means that this will likely be available in the future.

SpatialPixelsDataFrame

Given that a SpatialPixelsDataFrame has an attribute table we can use argument zcol = "column_name" to plot specific columns of the table. The zcol argument can be used with the …DataFrame versions of all Spatial* objects that are currently supported (see below).

mapView(meuse.grid, zcol = "dist")

Vector data objects

For spatial vector objects from the sp-package I would like to highlight the following arguments:

SpatialPoints(DataFrame)

data(meuse)
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")

# only one layer, all info in popups
mapView(meuse)

Another argument you can use with SpatialPointsDataFrame is radius which takes either an integer value (default is 10) or the character name of one of the columns in the attribute table. If the latter is provided, circles are scaled relative to the respective values in the column provided.

mapView(meuse, radius = "lead")

Apart from the radius argument, the same functionality is implemented for

SpatialPolygons(DataFrame)

data("gadmCHE")
mapView(gadmCHE)

and

SpatialLines(DataFrame)

data("atlStorms2005")
mapView(atlStorms2005, burst = TRUE)

Additional functionality

For the remote sensing inclined, mapview-package offers some additional functionality:

Red-Green-Blue images

viewRGB() provides the possibility to view Red-Green-Blue true-/false-color images of RasterStacks/Bricks. Any band combination is possible.

True color example

data("poppendorf")
viewRGB(poppendorf, 4, 3, 2) #Bands 4,3,2 correspond to Red Green Blue

Note, by default the color to render NA values is set to "transparent", so let’s change that to "black"

False color example

viewRGB(poppendorf, 5, 4, 3, na.color = "black") #Bands 5,4,3 correspond to NIR Red Green

Fun fact: within this tiny area one can find at least five breweries (see below)… Welcome to the land of Beer!

View extent only

As mentioned above, it is possible to view only the extent of Raster and bbox of Spatial objects. Clicking on the rectangle triggers a pop-up window showing the coordinates of the corners of the extent in epsg:4326.

Raster* object extent

viewExtent(poppendorf)

Spatial* object bbox

viewExtent(meuse, color = "black")

Overlays

To make mapview even more user-friendly, I have defined a +-method for combinations of

  • mapview + mapview
  • mapview + ANY (ANY here refers to all classes supported by mapView())
  • leaflet + ANY (so that one can easily pass spatial objects to existing leaflet maps)

This means you can easily combine objects by enabling things like

mapView(meuse.grid, zcol = "soil") + meuse

Basically, any combination is possible.

The zoom will be set to the extent of the objected that was added last. Also, layers are overlaid according to the order of calls. Note, however, that leaflet will always render vector data on top of raster data,

data("breweries91")
mapView(breweries91, color = "red") + viewRGB(poppendorf, 4, 3, 2) + viewExtent(poppendorf)

So here you have it, 5 breweries within the extent of the RGB image

leaflet integration

Integrating mapview and leaflet works both ways…

leaflet + mapview

You can first create a map using leaflet and then add spatial objects as you like with mapview allowing you to flexibly create maps in whatever way you like.

m <- leaflet() %>%
  addProviderTiles("Stamen.Toner")
m + meuse

mapview + leaflet

Given that objects created with mapView() contain a slot called map, you can first create a map using mapView() and then add components provided by leaflet. To reproduce the following example you will (at the time of this writing) need to install leaflet from the GitHub leaflet repository which provides new functionality to measure distances and areas. The following code will then add a control element in the top right corner of the map which lets you interactively create a polygon and will show its perimeter and area. Double-click to finish a measurement.

m <- mapView(breweries91) + viewExtent(poppendorf)
m@map %>% addMeasure(primaryLengthUnit = "meters", primaryAreaUnit = "sqmeters", 
                     activeColor = "#3D535D", completedColor = "#7D4479")

So now you can measure the distances between the breweries or figure out how big the area is that is covered by the extent…

Advanced example

As a final, more advanced, example of what can be done with mapview we will mimic a common workflow when working with spatial data. We will calculate the mean January temperature for each district of Switzerland using raster::extract. Then we will render this on top of the mean January temperatures and also add a hill-shade raster of Switzerland. All data sets are included in mapview. They were created following the examples shown in the help pages of raster::getData and raster::hillShade. This example also highlights a few arguments that can be adjusted which were not mentioned so far.

data("tmeanCHE")
data("hillshadeCHE")
data("gadmCHE")

t <- extract(tmeanCHE, gadmCHE, fun = mean, na.rm = TRUE, sp = TRUE)
t$Jan.mean.T <- round(t$Jan.mean.T, 2)

m1 <- mapView(hillshadeCHE, color = grey.colors(3), 
              layer.opacity = 1, layer.name = "Hillshade",
              map.types = "Acetate.hillshading", legend = FALSE)
m2 <- mapView(tmeanCHE, layer.opacity = 0.6, 
              map.types = c("OpenStreetMap.BlackAndWhite", 
                            "OpenStreetMap.HOT"),
              layer.name = "JanTemp")
m3 <- mapView(t, zcol = "Jan.mean.T", fillOpacity = 0.9)
m1 + m2 + m3

Note how we

I hope this highlights that despite providing sensible default behaviour for viewing spatial data mapview also provides a fair deal of flexibility.

And for those who…

… still wonder where exactly the difference between mapview and leaflet lies, consider the following:

Let’s reproduce a simple mapView() call using leaflet()

Here’s the mapview version:

mapView(meuse)

And now let’s get the identical result using leaflet:

## first we need to reproject meuse to geographic coordinates
meuse <- spTransform(meuse, 
                     CRSobj = CRS("+proj=longlat +datum=WGS84 +no_defs"))

## then we need to generate the text for the pop-ups
df <- as.data.frame(sapply(meuse@data, as.character),
                    stringsAsFactors = FALSE)

nms <- names(df)

txt_x <- paste0("x: ", round(coordinates(meuse)[, 1], 2))
txt_y <- paste0("y: ", round(coordinates(meuse)[, 2], 2))

txt <- rbind(sapply(seq(nrow(meuse@data)), function(i) {
  paste(nms, df[i, ], sep = ": ")
}), txt_x, txt_y)

txt <- sapply(seq(ncol(txt)), function(j) {
  paste(txt[, j], collapse = " <br/> ")
})

## finally we can create our map
leaflet(meuse) %>% 
  addProviderTiles("OpenStreetMap", group = "OpenStreetMap") %>% 
  addProviderTiles("Esri.WorldImagery", group = "Esri.WorldImagery") %>%
  addCircleMarkers(lng = coordinates(meuse)[, 1],
                   lat = coordinates(meuse)[, 2],
                   group = "meuse",
                   color = mapViewPalette(3)[length(mapViewPalette(3))],
                   popup = txt) %>%
  addLayersControl(position = "bottomleft",
                   baseGroups = c("OpenStreetMap",
                                 "Esri.WorldImagery"),
                   overlayGroups = "meuse")

I guess this highlights quite well what mapview is intended for. It serves to save you some effort to interactively examine your spatial data by giving you the possibility to render it using a one-liner of code.

Final thoughts

Where to go from here?

Well, there are still a lot of limitations that need addressing:

In future releases I would like to

I hope that mapview will prove useful for some people out there.

If you have any feedback, please don’t hesitate to contact me.

Bug reports should be filed at https://github.com/environmentalinformatics-marburg/mapview/issues

Best,

Tim